function [X1_control,Y1_control,X2_control,Y2_control] = plot_2_BForm_2(flag,dist,cin1,cin2,v1,v2,v3,v4,d1,d2)
% [X1_control,Y1_control,X2_control,Y2_control] =
% plot_2_BForm(flag,dist,cin1,cin2,v1,v2,v3,v4,d) 
% which return the coordinates of control points for two surfaces
if nargin <= 0
    flag = 2;
end
if nargin <= 1
    dist = 0.15;
end
if nargin <= 2
    d1 = 3;
    d2 = 4;
    v1 = [   2 , 1.5   ];
    v2 = [   0 , 1.732 ];
    v3 = [   1 , 0     ];
    v4 = [-1.5 , 0     ];
    cin1 = [3 , 4 4.2 , 4.5 4.8 4.5 , 3.9 4.6 4.8  4.4   ]';
    cin2 = [3 , 4 4.2 , 4.4 4.8 4.5, 4.6 5 5.1 4.4 , 3.9 4.4 4.8 4.5 4.0 ]';
end

%T1 = [v1,v2,v3];
%T2 = [v4,v3,v2];

%desc = desc_pattern(d1+d2); % use not so constant ,so generate here is not so costy
[cr,desc,asce] = get_auxillary_mat(d1+d2);
plot_d = 20; flag=3;

% calculate the data for surface
[I,J,K] = indices(plot_d);   % for sample point calculating
I = I/plot_d; J = J/plot_d; K = K/plot_d;
[phi1,dphidxi1,dphideta1] = eval_bb(d1,I',J'); % get the basis function value on sample point
[phi2,dphidxi2,dphideta2] = eval_bb(d2,I',J'); % get the basis function value on sample point
% until here, phi is sized m1 x m2
% m1 = (d+1)*(d+2)/2; m2 = (plot_d+1)*(plot_d+2)/2

tri = make_mesh(plot_d);  % for function trimesh, it it common for all triangles

% x-y coordinates for triangle T1 and T2
X1 = I*v1(1) + J*v2(1) + K*v3(1);
Y1 = I*v1(2) + J*v2(2) + K*v3(2);
% dist = 0.15;
X2 = I*v4(1) + J*v3(1) + K*v2(1) - dist;
Y2 = I*v4(2) + J*v3(2) + K*v2(2) - dist;

Z1 = phi1' * cin1 ; % size(Z) = [m2,1]
Z2 = phi2' * cin2 ; % size(Z) = [m2,1]

% calculate the data for control polygon and triangle
tri_control_1 = make_mesh(d1);  % for function trimesh, it it common for all triangles
[I_control_1,J_control_1,K_control_1] = indices(d1);   % for sample point calculating
I_control_1 = I_control_1/d1; J_control_1 = J_control_1/d1; K_control_1 = K_control_1/d1;
tri_control_2 = make_mesh(d2);  % for function trimesh, it it common for all triangles
[I_control_2,J_control_2,K_control_2] = indices(d2);   % for sample point calculating
I_control_2 = I_control_2/d2; J_control_2 = J_control_2/d2; K_control_2 = K_control_2/d2;

% x-y coordinates for triangle T1 and T2
X1_control = I_control_1*v1(1) + J_control_1*v2(1) + K_control_1*v3(1);
Y1_control = I_control_1*v1(2) + J_control_1*v2(2) + K_control_1*v3(2);
X2_control = I_control_2*v4(1) + J_control_2*v3(1) + K_control_2*v2(1) - dist;
Y2_control = I_control_2*v4(2) + J_control_2*v3(2) + K_control_2*v2(2) - dist;

% begin to draw 
figure; set(gcf,'Position',[100 100 600 400]); hold on;
if flag == 1 | flag == 3 % then draw the surface
    trisurf(tri,X1,Y1,Z1);    % draw surface 1
    trisurf(tri,X2,Y2,Z2); 
    shading interp;  % draw the surface
end
if flag == 2 | flag == 3  % then draw the control polygon 
    trisurf(tri_control_1,X1_control,Y1_control,cin1);                % draw the control polygon
    trisurf(tri_control_2,X2_control,Y2_control,cin2);                % draw the control polygon
end
% at last the grid
trisurf(tri_control_1,X1_control,Y1_control,zeros(size(X1_control)));
trisurf(tri_control_2,X2_control,Y2_control,zeros(size(X2_control)));

axis([-1.7 2 -0.5 1.8 0 6]);
view(pi/4,30); hold off; % axis off;